#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <cctype>

// —— 安全 malloc —— //
static void* xmalloc(size_t n) {
    void *p = std::malloc(n);
    if (!p) {
        std::fprintf(stderr, "内存分配失败（%zu 字节）。数据维度可能过大。\n", n);
        std::exit(2);
    }
    return p;
}

static int is_line_blank(const char *s) {
    while (*s) { if (!std::isspace((unsigned char)*s)) return 0; ++s; }
    return 1;
}

// —— 两遍读取 CSV：先数行列，再读入 —— //
static int read_csv_two_pass(const char *filename, double **X_out, int *rows_out, int *cols_out) {
    std::FILE *fp = std::fopen(filename, "r");
    if (!fp) { std::perror("打开CSV失败"); return -1; }

    char *line = (char*)xmalloc(1<<16);
    int rows = 0, cols = -1;

    // 第一遍：统计行列
    while (std::fgets(line, 1<<16, fp)) {
        if (is_line_blank(line)) continue;
        int c = 0;
        for (char *tok = std::strtok(line, ",\r\n"); tok; tok = std::strtok(nullptr, ",\r\n")) c++;
        if (c == 0) continue;
        if (cols == -1) cols = c;
        else if (c != cols) {
            std::fprintf(stderr, "第 %d 行列数 %d 与首行 %d 不一致。\n", rows+1, c, cols);
            std::fclose(fp); std::free(line); return -2;
        }
        rows++;
    }
    if (rows == 0 || cols <= 0) { std::fclose(fp); std::free(line); return -3; }

    if (std::fseek(fp, 0, SEEK_SET) != 0) { std::fclose(fp); std::free(line); return -4; }

    double *X = (double*)xmalloc(sizeof(double) * rows * cols);
    int r = 0;
    while (std::fgets(line, 1<<16, fp)) {
        if (is_line_blank(line)) continue;
        int c = 0;
        for (char *tok = std::strtok(line, ",\r\n"); tok; tok = std::strtok(nullptr, ",\r\n")) {
            X[r*cols + c] = std::atof(tok);
            c++;
        }
        if (c > 0) r++;
    }
    std::fclose(fp); std::free(line);

    *X_out = X; *rows_out = rows; *cols_out = cols;
    return 0;
}

// —— 预处理 —— //
static void center_columns(double *X, int rows, int cols) {
    for (int j = 0; j < cols; ++j) {
        double s = 0.0;
        for (int i = 0; i < rows; ++i) s += X[i*cols + j];
        double mean = s / rows;
        for (int i = 0; i < rows; ++i) X[i*cols + j] -= mean;
    }
}
static void standardize_columns(double *X, int rows, int cols) {
    for (int j = 0; j < cols; ++j) {
        double s = 0.0, sq = 0.0;
        for (int i = 0; i < rows; ++i) {
            double v = X[i*cols + j];
            s += v; sq += v*v;
        }
        double mean = s / rows;
        double var  = std::fmax(1e-12, sq / rows - mean*mean);
        double stdv = std::sqrt(var);
        for (int i = 0; i < rows; ++i) X[i*cols + j] = (X[i*cols + j] - mean) / stdv;
    }
}

// —— 协方差：cov = (1/(n-1)) X^T X —— //
static void compute_covariance(const double *X, double *cov, int rows, int cols) {
    double scale = 1.0 / (rows - 1);
    for (int a = 0; a < cols; ++a) {
        for (int b = a; b < cols; ++b) {
            double s = 0.0;
            for (int i = 0; i < rows; ++i) s += X[i*cols + a] * X[i*cols + b];
            cov[a*cols + b] = cov[b*cols + a] = s * scale;
        }
    }
}

// —— Jacobi 特征分解（堆分配，不用大栈） —— //
static void eye(double *V, int n) {
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            V[i*n + j] = (i == j) ? 1.0 : 0.0;
}
static void jacobi_eigen(double *A, double *V, double *d, int n, int max_iter, double tol) {
    eye(V, n);
    for (int iter = 0; iter < max_iter; ++iter) {
        int p = 0, q = 1;
        double maxv = std::fabs(A[p*n + q]);
        for (int i = 0; i < n; ++i)
            for (int j = i + 1; j < n; ++j) {
                double v = std::fabs(A[i*n + j]);
                if (v > maxv) { maxv = v; p = i; q = j; }
            }
        if (maxv < tol) break;

        double app = A[p*n + p], aqq = A[q*n + q], apq = A[p*n + q];
        double phi = 0.5 * std::atan2(2.0 * apq, (aqq - app));
        double c = std::cos(phi), s = std::sin(phi);

        for (int k = 0; k < n; ++k) {
            double aik = A[p*n + k], aqk = A[q*n + k];
            A[p*n + k] = c*aik - s*aqk;
            A[q*n + k] = s*aik + c*aqk;
        }
        for (int k = 0; k < n; ++k) {
            A[k*n + p] = A[p*n + k];
            A[k*n + q] = A[q*n + k];
        }
        A[p*n + p] = c*c*app - 2*s*c*apq + s*s*aqq;
        A[q*n + q] = s*s*app + 2*s*c*apq + c*c*aqq;
        A[p*n + q] = A[q*n + p] = 0.0;

        for (int k = 0; k < n; ++k) {
            double vip = V[k*n + p], viq = V[k*n + q];
            V[k*n + p] = c*vip - s*viq;
            V[k*n + q] = s*vip + c*viq;
        }
    }
    for (int i = 0; i < n; ++i) d[i] = A[i*n + i];
}

// —— 排序（同步重排特征向量列） —— //
static void sort_eigenpairs_desc(double *d, double *V, int n) {
    for (int i = 0; i < n - 1; ++i) {
        int m = i;
        for (int j = i + 1; j < n; ++j)
            if (d[j] > d[m]) m = j;
        if (m != i) {
            std::swap(d[i], d[m]);
            for (int r = 0; r < n; ++r)
                std::swap(V[r*n + i], V[r*n + m]);
        }
    }
}

// —— 矩阵乘法 —— //
static void matmul(const double *A, const double *B, double *C, int m, int n, int k) {
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < k; ++j) {
            double s = 0.0;
            for (int t = 0; t < n; ++t) s += A[i*n + t] * B[t*k + j];
            C[i*k + j] = s;
        }
    }
}

// —— 写 CSV —— //
static int write_matrix_csv(const char *filename, const double *M, int rows, int cols) {
    std::FILE *fp = std::fopen(filename, "w");
    if (!fp) { std::perror("写文件失败"); return -1; }
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            std::fprintf(fp, "%.10g", M[i*cols + j]);
            if (j + 1 < cols) std::fputc(',', fp);
        }
        std::fputc('\n', fp);
    }
    std::fclose(fp);
    return 0;
}
static int write_vector_csv(const char *filename, const double *v, int n) {
    std::FILE *fp = std::fopen(filename, "w");
    if (!fp) { std::perror("写文件失败"); return -1; }
    for (int i = 0; i < n; ++i) std::fprintf(fp, "%.10g\n", v[i]);
    std::fclose(fp);
    return 0;
}

// —— 主程序 —— //
int main(int argc, char **argv) {
#ifdef _WIN32
    system("chcp 65001 > nul");
#endif

    const char *input = "data.csv";
    int K = 2;
    int use_standardize = 1;

    if (argc >= 2) input = argv[1];
    if (argc >= 3) K = std::atoi(argv[2]);
    if (argc >= 4) use_standardize = std::atoi(argv[3]);

    std::printf("读取：%s ...\n", input);
    double *X = nullptr; int rows = 0, cols = 0;
    if (read_csv_two_pass(input, &X, &rows, &cols) != 0) {
        std::fprintf(stderr, "读取失败：请确认文件存在且无表头、逗号分隔。\n");
        return 1;
    }
    if (rows < 2 || cols < 1) { std::fprintf(stderr, "数据规模不足。\n"); std::free(X); return 1; }
    if (K <= 0 || K > cols) { K = (cols >= 2) ? 2 : 1; }
    std::printf("数据规模：%d × %d，目标维度 K=%d\n", rows, cols, K);

    std::printf("%s...\n", use_standardize ? "标准化" : "中心化");
    if (use_standardize) standardize_columns(X, rows, cols);
    else                 center_columns(X, rows, cols);

    // 注意：协方差是 cols×cols，如果列非常大（比如 > 3000），内存会很大
    size_t cov_bytes = sizeof(double) * (size_t)cols * (size_t)cols;
    std::printf("分配协方差矩阵内存：%.2f MB\n", cov_bytes / (1024.0*1024.0));
    double *cov = (double*)xmalloc(cov_bytes);

    std::printf("计算协方差...\n");
    compute_covariance(X, cov, rows, cols);

    double *V    = (double*)xmalloc(sizeof(double) * (size_t)cols * (size_t)cols);
    double *eigs = (double*)xmalloc(sizeof(double) * (size_t)cols);

    std::printf("Jacobi 分解...\n");
    jacobi_eigen(cov, V, eigs, cols, 500, 1e-9);

    std::printf("排序...\n");
    sort_eigenpairs_desc(eigs, V, cols);

    double sum_all = 0.0, sum_k = 0.0;
    for (int i = 0; i < cols; ++i) if (eigs[i] > 0) sum_all += eigs[i];
    for (int i = 0; i < K;    ++i) if (eigs[i] > 0) sum_k   += eigs[i];
    std::printf("信息保留率（前 %d 个主成分）：%.2f%%\n", K, sum_all>0? (sum_k/sum_all*100.0):0.0);

    std::printf("计算投影...\n");
    double *W = (double*)xmalloc(sizeof(double) * (size_t)cols * (size_t)K);
    for (int j = 0; j < K; ++j)
        for (int i = 0; i < cols; ++i)
            W[i*K + j] = V[i*cols + j];

    double *Z = (double*)xmalloc(sizeof(double) * (size_t)rows * (size_t)K);
    matmul(X, W, Z, rows, cols, K);

    write_matrix_csv("reduced.csv", Z, rows, K);
    write_vector_csv("eigenvalues.csv", eigs, cols);
    write_matrix_csv("eigenvectors.csv", V, cols, cols);

    std::printf("完成：reduced.csv / eigenvalues.csv / eigenvectors.csv\n");

    std::free(X); std::free(cov); std::free(V); std::free(eigs); std::free(W); std::free(Z);
    return 0;
}